Magnetic droplets in a metal close to a ferromagnetic quantum critical point 
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Using analytical and path integral Monte Carlo methods, we study the susceptibility Xdc{T) of 
a spin-S impurity with XY rotational symmetry embedded in a metal. Close to a ferromagnetic 
quantum critical point, the impurity polarizes conduction electrons in its vicinity and forms a large 
magnetic droplet with moment M ^> S. At not too low temperatures, the strongly damping param- 
agnon modes of the conduction electrons suppress large quantum fluctuations (or spin flips) of this 
droplet. We show that the susceptibility follows the law xM T ) = (M 2 /T)[l - (irg) -1 ln(gE /T)], 
where the parameter j > 1 describes the strong damping by conduction electrons, and Eo is the 
bandwidth of paramagnon modes. At exponentially low temperatures TcT, ~ Eo exp(— 7rg/2) we 
show that spin flips cannot be ignored. In this regime we find that Xdc(T) fa Xde(0)[l— (2/3)(T/T*) 2 ], 
where Xdc(0) ~ M 2 /T* is finite and exponentially large in g. We also discuss these effects in the 
context of the multi-channel Kondo impurity model. 
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I. INTRODUCTION 



Magnetic impurities in a nearly ferromagnetic Landau 
Fermi-liquid can induce large magnetic droplets by po- 
larizing conduction electrons in their vicinityi The large 
magnetic polarizability of the conduction electrons can 
be described in terms of low-energy collective excitations, 
paramagnons^ An impurity spin dressed by these soft 
paramagnon excitations forms a magnetic droplet, and 
the size of such droplets, determined by the spatial dis- 
persion of the paramagnons, can greatly exceed typical 
interatomic distances in the proximity of the quantum 
critical point. The dynamics of such a droplet are essen- 
tially determined by the paramagnon modes which damp 
the orientational motion of the droplet. At not too low 
temperatures, the fluctuations of the droplet's moment 
are small if the damping of the angular motion is strong. 
As the temperature is lowered, due to these small fluc- 
tuations the effective damping decreases slowly, usually 
as a power law or logarithmically (see below). This de- 
crease of the effective damping is reflected in the temper- 
ature dependence of the impurity's magnetic susceptibil- 
ity which increases at a rate slower than the Curie- Weiss 
law, x<ic{T) ~ T , as the temperature is decreased. 

The relevance of large quantum fluctuations (or spin 
flips) of these overdamped magnetic droplets, the main 
interest of this paper, has been a topic of active inves- 
tigation recently. Consensus in this matter has proved 
elusive. The various existing points of view seem to 
agree that the susceptibility of large overdamped droplets 
at not very low temperatures should obey Xdc{T) ~ 
T~ 1+Q , where a < I is non-universal. Millis, Morr and 
Schmalian 3 found that in the case of a magnetic defect 
with Ising symmetry, quantum tunneling is suppressed 
altogether, and that the power-law temperature depen- 
dence of XdciT) extends down to T = K. Furthermore, 
they indicated that their conclusion was relevant even for 
defects with a continuous symmetry 4 Castro-Neto and 



Jones^, in an earlier work, analysed the same problem 
with both ferromagnetic and antiferromagnetic clusters. 
In a broader context, our analysis forms a part of the 
general problem of disorder in nearly quantum critical 
metals. Three decades ago, Griffiths^ and McCoy?£ pre- 
dicted a non-analytic temperature dependence of mag- 
netization in nearly-ordered ferromagnetic Ising models 
with bond disorder. Throughout the years after that, the 
research on quantum Griffiths, Kondo disorder and local 
criticality 8 '^ problems remains of wide interest. 

We study the magnetic susceptibility, Xdc(T), of a spin- 
s' impurity with XY rotational symmetry coupled to the 
conduction electrons through an exchange interaction, 



H ex = JS- ct(0)<rc(0), 



(1) 



in a metal close to a ferromagnetic quantum critical 
point. We employ both analytical and path integral 
Monte Carlo (PIMC) techniques. The two main results 
of our analysis are as follows. First, we show that at not 
very low temperatures, due to small quantum fluctua- 
tions of the XY droplet the susceptibility evolves as 



Xdc(T) 



■ng V T 



We shall show that the result is valid over an expo- 
nentially large range of temperature, gEo 3> T 3> 
Eq exp(— ng/2). Here the parameter g fa ( Jn(e_p)) 2 /(f + 
F a ) ^> I represents the damping of the droplet's fluctu- 
ations, and \M\ w \J\n(e F )S/(l + F a ) > S is the mag- 
netic moment of the droplet. F a — > —I is the Landau 
Fermi-liquid parameter denoting distance from the quan- 
tum critical point. To leading order in g^ 1 , the logarith- 
mic form and the power law, Xdc{T) = (r/T )" 1+1/(7r9) , 
are the same. However at lower temperatures, the two 
expressions can differ significantly. Second, and more 
significantly, we show that at exponentially low temper- 
atures T « T, = Eq exp(— ng/2), full 2tt spin rotations 



2 



cannot be ignored. The susceptibility at T = saturates 
in this model, and at finite temperatures T<T„ 



Xdc(T) w Xdc(0)[l — (2/3)(T/T*) 



where Xdc(O) ~ M 2 /T* is of the order of the susceptibility 
at T ~ T*. The zero-temperature susceptibility Xdc(O) is 
exponentially large in g. The conclusion is that for the 
XY magnetic defect, spin flips are important at low tem- 
peratures, and their effect is to remove the divergence of 
susceptibility at T = K. The low-temperature behav- 
ior we obtain for the impurity susceptibility is similar to 
that seen in the usual Kondo problemifi*ii for J > 0. 
The regime of proliferation of spin flips below T < T* 
is analogous to the strong-coupling regime of the Kondo 
problem below the Kondo temperature. In contrast to 
the usual Kondo effect where the sign of the exchange 
coupling between the impurity spin and conduction elec- 
trons, </, is important, our results are independent of the 
sign of J. We consider here only the effects dependent on 
the damping g and therefore on even powers of J. In the 
closing section we discuss further the experimental real- 
izations and differences of magnetic droplet phenomena 
from the standard Kondo effec1iii*i£. 

We assume that the coupling between a spin and elec- 
trons in a single channel is small, Jn(e F ) <C 1. There- 
fore the terms beyond the Born approximation can be 
neglectedi Importantly, the overall coupling constant 
g « (Jn(e F )) 2 / (1 + F a ) is large in the proximity of the 
quantum critical point F a — > — 1 . The magnetic moment 
of the droplet, M = (1 - Jn(e F )/(l + F a ))S, is large, 
because the droplet is dressed by (or coupled to) a large 
number of electron channels, N c h = 1/(1 + F a ). Close 
to the critical point, the contribution of such magnetic 
droplets to the susceptibility and resistivity can over- 
shadow other impurity effectsi 

The weak damping (g <C 1) limit of the same problem 
was first studied by Larkin and Melnikov (LM)i. Both 
resistivity and susceptibility were found to decrease log- 
arithmically with temperature irrespective of the sign of 
exchange coupling J as opposed to the usual Kondo effect 
where resistivity increases logarithmically with tempera- 
ture if J > 0. 

The rest of this paper is organized as follows. In Sec ITU 
we derive our model for dissipative dynamics of the mag- 
netic droplet beginning with the exchange Hamiltonian, 
Eq.Q. The correlation function and susceptibility of 
the droplet are studied analytically in Sec lIIII We study 
both not-too-low temperatures where tunneling effects 
are negligible, and very low temperatures where tunnel- 
ing makes important contributions. Sec lIVI provides de- 
tails of the numerical path integral Monte Carlo (PIMC) 
method we employ. In Sec[V] we present the results of 
our analytic and numerical study. Scc lVIl contains a dis- 



II. MODEL AND FORMALISM 

The magnetic properties of the metal are determined 
by the spin susceptibility of the conduction electrons, 



£>(k,w„)<5« 



dre 



ik-r 







dre l "" T (CT l (0,0)(T J (r,T)), 



where a l (r,r) is the conduction electron spin density, 

^(r,r) =cj 4 (r,r)cr* w c 1 ,(r,r), 

and u> n = 2nTn. The static part of the susceptibility is 
D(0, 0) = 2n(e F ) /(1 + F a ) in terms of the standard Lan- 
dau Fermi- liquid parameter F a . Close to a ferromagnetic 
instability, F a ss —I, due to the large static magnetic sus- 
ceptibility the impurity induces a large magnetic droplet 
with effective moment M = S[l — Jn(e F )/ (1+F a )]. In the 
case of an antiferromagnetic exchange coupling (J > 0), 
the droplet 's magnetic moment is polarized in the op- 
posite direction to the impurity spin, while in the case 
of ferromagnetic coupling (J < 0) the droplet's magnetic 
moment and the impurity spin are locked parallel to each 
other. The low- lying magnetic excitations of the conduc- 
tion electrons (which also constitute the major part of 
the droplet) are strongly damped, as can be seen in the 
expression for spin susceptibility at small k and u> n , 



D(k,uj r 



2n{e F ) 



(&k) : 



7T \u>n\ 
2 k v F 



(2) 



£o is a length scale of the order of interatomic distances. 
The dynamics of the impurity are determined by the local 
susceptibility D{uj n )r^ 

D(cj n )-D(0) = J(d 3 k) [D(k,u n )-D(k,0)] 
. Trn(e F ) 2 \(v n \ 

~ 2{s, Q k F f(i + F a y y) 

Eq.(J3J) is valid at low enough frequencies, 

\u n \ «£ = eF(l + F Q ) 3 / 2 . 

At higher frequencies, E <g; \uj n \ <g; e F , Eq.© should be 
replaced with 



33/2(£ fc F )8/3 



e F 



Integrating out the conduction electrons (in the per- 
turbation series of (Jn(e F )) 2 <C I) in Eq.Q) results in a 
dissipative action for the impurity^ 



j 2 

[S] = — / drdr'Sir) ■ S(t')D(t - t'). 
2 Jo 



Simp fSl 



We are interested in the impurity dynamics at low tem- 
perature, so using Eq.JSJ) for the interaction, 

_ . . irgT 2 [P , , 1 - n(r) • n(r') 
S imp [n] = -2— / drdr' 2l _J,_ _ )N V (5) 



n 



sin 2 (7rT(r-r')) 
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In Eq.©, S(r) = Sn(r), n(r) 2 
the damping g given by 



1, and we assume that 



9 



2(1 + F a ) 



[JSn{e F )f 



(6) 



is large (g > 1). This is always possible sufficiently close 
to the transition. The form of the interaction we chose in 
Eq.(|5j) is valid only up to an energy Eq. We may impose 
this cutoff through an additional regularizing term in the 
action, 



1 f p 



(7) 



This method of imposing the cutoff is not unique. For 
instance, an equally valid option would have been to 
introduce a short time cutoff, r c in the interaction, 
gT 2 / sin 2 (ttTt) -> gT 2 / sin 2 (7rT^r 2 + r c 2 ). The cutoff 
appears in the correlation functions only a function of 
Eqt or t/t c , and as long as gE$T 3> 1, the results are in- 
dependent of the manner in which the cutoff is imposed. 
Parametrizing n(r) = (cos tp{r), sin (p(r)), our dissipative 
model for the XY impurity takes the final form, 



«%] = t^t / dr [d r (p) 

4-C/o Jq 



TrgT 2 



dr dr 



M / drh_L(r) • n(r) + 
Jo 

,1 - cos(<p(t) - (p(t')) 



sin 2 (7rT(r - r')) 



(8) 



where M = S[l — Jti(cf) / {1+ F a )] is the bare droplet mo- 
ment, and hj^ is an in-plane magnetic field. The Matsub- 
ara fields ip(r) satisfy periodic boundary conditions up to 
multiples of 2ir: ip(r + 0) = (p(r) + 2irk, where k is an 
integer called the winding number. The full phase <^(t) 
can always be written in the form tp(r) = 2irTkr + </>(r), 
where the residual phase obeys periodic boundary condi- 
tions 4>(t + [3) — (p(r). 

We are mostly concerned here with the rotational (or 
orientational) motion of the droplet 's moment. We as- 
sume throughout the main text (but see the concluding 
section) that the exchange coupling is sufficiently strong 
to suppress the fluctuations of the amplitude of the bare 
magnetic moment M , namely, in the quantum limit or at 
low temperatures, JS » kT as well as E ^S> kT. It also 
suits us to regularize the tunneling action through the in- 
troduction of a kinetic term because the model in Eq. JSJ) 
appears in numerous contexts. We have already men- 
tioned the recent work by Millis and co-workers^ in which 
the authors studied the dynamics of magnetic defects in 
nearly quantum-critical metals, where the defects are re- 
gions of ordered phase formed due to a local enhance- 
ment of the transition temperature. For defects with XY 
symmetry, they arrived at essentially the same strongly- 
damped model as Eq. iJSJ , although in their model the ki- 
netic term did not appear simply as a means for imposing 
a cutoff but had a definite physical meaning as the contri- 
bution to the droplet action from the magnon part of the 
dispersion curve. Our analysis should be valid for such 



systems as well. The action in Eq.© also arises from 
an Ambegaokar-Eckern-Schoni^ treatment of tunneling 
through a quantum dot. There, the physical meaning of 
Eq is the charging energy for the quantum dot. 

In this paper we study the impurity spin correlator 
C(r), 

C(r) = (co S (<p(r)-<p(0))), (9) 
and the zero- frequency impurity susceptibility Xdc{T), 

Xdc(T) = M 2 / dr (n(r) • n(0)> 



M 



dr (cob(<p(t) - <p(0))). (10) 



Calculation of the imaginary part of the susceptibility as 
well as transport properties like resistivity involve sub- 
tleties associated with analytic continuation to real fre- 
quencies. These will be studied in a later work. In this 
paper we calculate only the real part of the susceptibility 
as shown in Ea. ((TU)l . 



III. ANALYSIS OF IMPURITY SPIN 
CORRELATION FUNCTION AND 
SUSCEPTIBILITY, g > 1 

A large value of g tends to suppress large fluctua- 
tions (tunneling) of the droplet moment. Physically, the 
droplet couples to a large number of channels {g oc N c h) 
of the conduction-electron continuum, which makes spin 
flips difficult. We show below that this is not the case 
at very low temperatures, where spin flips (we use the 
terminology of spin flips and tunneling interchangeably 
in order to discuss XY and Ising symmetry simultane- 
ously) can occur even for large g. In the first part of this 
section, we consider not-too-low temperatures where tun- 
neling may be disregarded. In the latter half of this sec- 
tion, we discuss the limits beyond which tunneling may 
not be ignored, and analyze the effect of tunneling on the 
correlation function and susceptibility. 



A. Ignoring winding numbers 

Let us begin by studying the action, Eq. © (with 
hj^ = 0), by ignoring winding numbers (k — 0) and ex- 
panding cos(</j(t) — ^>{t')) to quartic order in the phase 
difference. We then gather the Gaussian terms in the 
resulting action as the bare term and treat the quartic 
term as an "interaction" (S = S sauss + S iDt ) : 

4£o Jo 
ngT 2 rP 



4 

■KgT 
48 



dT{d T <t>) 2 + 

, W>(t)-0(t')) 2 



drdr' 



sin 2 (7i-T(r-T'))' 



20 '"drdr' W-)-^)) 4 



o 



sin^(7rT(r- t')) 



(11) 



(12) 
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The bare term may be diagonalized by going over to the 
frequency representation, 



Essentially, fbareM is logarithmic in r, with a lower cut- 

o n r\ an in 

gE 



off and an upper cutoff 



(r) = J^(a„ cosw„r + 6„sinu; n r), thus, 



n=l 

oo 



Tvr 2 



Er 



-n + irgn 



In the following discussion, we will need the bare Green 
function, which is given exactly by the following sum: 



^barc(r) = ((0(r) - 0(O)) 2 ) barc 
_ „ 1 - coso; k t 

^! ^r" + van 



(13) 



-Fbarc(T) 



2^o Irl 



r < 



9^0 



2_ i 2e^ g £o sin 7tTt 
■Kg 7rT 5 



T > 



1 



(14) 



By 



Consider now the renormalization of S sanss by S lnt 
contracting two of the four fields appearing in S sanss us- 
ing the bare Green function, we obtain a one-loop renor- 
malization of the bare action by the interaction: 



irgT 



2 r P 



drdr 



, WT)-^(r')) 2 



jo sin 2 (7rT(r-r')) 
ivgT 2 — ^ 



((0(r) - ^(r')) 2 )b arc 



/" 1 s 

E / drdr ' — ^ 



,l-cosK(r-r')) v (0(r) - 0(r')) 2 

X 



'-n 2 +irgn sin (-kT{t — r')) 



The effective Gaussian action, S" c S ff USS [</>] = S' gauss [0] + 
(5S ,gauss [^)], with fluctuations considered up to one loop is 
therefore 



q. 



off M^EK 

71=1 



7T 2 T 



where the effective coupling to one loop, 

gven(n) =g[l ln ^T 

V irg nl 



(15) 



(16) 



is smaller than the bare coupling g, and the reduction 
is strongest at low n. Since g oc M, this means that 
low-energy paramagnon fluctuations of the conduction 
electrons that constitute a large fraction of the droplet 
cause a logarithmic reduction of its moment. 

The phase correlation function with the dressed cou- 
pling g rcn (n), 

C(r) = (cos(0(r) - 0(0))) = exp(-((0(r) - 0(O)) 2 )/2) 

1 — cos(w n r) 



exp 



E 



ifnH^n l-iln^ 



nT 



1 



1 _ JL ln( 5 S Q |r|) + 0(.g- 3 ln 3 ( 5j B |r|)), (17) 
irg 



differs at 0(g 2 ) from the bare correlator 



Cbarc(T) = exp(-F ba rc(T)/2) ~ (.9^0 M) . 



(18) 



At large |r|, the difference between the bare correlator 
and dressed correlator can be substantial. Eq. l|17l) is a 
one-loop calculation for C(r). Notice that if the effective 
g ren has only a logarithmic correction as in Ea- ljlfijl . the 
correlator C(r) has the same logarithmic correction as in 
Ea. (|17|) . To obtain the next 0(g~ 1 ) contribution to the 
coupling g ren in the Gaussian approximation, we need 
to calculate the two-loop diagrams. These include only 
the second-order diagrams from the quartic interaction 
in Ea. (|12fl . but also one diagram from the sixth-order 
fluctuation expansion of the tunneling term in Eq.©. 
One may check that the most singular contribution in 
both classes of diagrams is 0{g~ 1 In 2 (gE/T)). What 
is also clear that all other diagrams are of the order 
0(g 2 ~ m ~ l \n m (gE/T)), where m, I are natural numbers. 
Therefore if ln(gE/T) <C g, Eqs. (|16ll7H are good approx- 
imations. Since the logarithm is a slow function of r or 
T, there is a large range of r or T where C(r) can be 
approximated by Ea. (|17|) . Such a perturbative analy- 
sis does not preclude a qualitatively different behavior of 
C(t) at long r's, which indeed occurs as we discuss in the 



5 



next section. Ea. (|17fl leads to a magnetic susceptibility, 



Xdc{T) 



M 2 



iln^ 

irg T 



(19) 



In deriving Eq. (fTU|l . we assumed that C(r) is effectively 
constant over the whole of the interval [0,0], and that 
the errors in this assumption simply modify the cutoff of 
the logarithm. 

At very large |r| > 1/T gauss ~ {gE^e^, C(r) in 
Ea. ((TT|) becomes unphysically negative, which sets a low- 
temperature limit of validity of our perturbative analysis. 
Similarly, Eq. l|19|) is unable to resolve the question as to 
how the susceptibility behaves as T — > since our present 
perturbation theory, which ignores winding numbers and 
higher-order residual fluctuations, breaks down at tem- 
peratures T < T gauss = gE e^ 7rg . Near T = T gauss , 
Ea. (|19|l gives an exponentially large susceptibility, 



XdciTgauss) 



M 2 (e 7rs -1-wg) 
ng 2 E 



At lower temperatures T < T gauss one must consider 
winding numbers, higher-order residual fluctuations, and 



non-perturbative (in 1/g) contributions to Xdc- Our nu- 
merical calculations show that for g > 1, C(r) visibly 
falls below Eq.ljrT} around C(r) w 1/2, so Eq.fCJ be- 
comes inaccurate even before |r| « l/T gauss is reached. 

We emphasize that the logarithmic temperature de- 
pendence in Eo. (|19|l is not an approximate expansion of 
^(T/gEo) 1 /^ 9 ^ to order 1/g, but is in fact more accu- 
rate than this, in the range gEo > T> T gauss . 



B. Taking winding numbers into account 

We now turn our attention to paths fir) = UkT + ififr) 
with a finite winding number fc. We have chosen to ex- 
pand the full phase ip in terms of residual-phase fluctua- 
tions, (j>, about the 'classical' paths <p c1 (t) = lo^t. There 
are reasons for doing so. Firstly, these classical paths are 
solutions of the Euler-Lagrange equations, so they are 
stationary points of the action Sk [4>] ■ Secondly, expand- 
ing the action Eq.JHJ to second order in the fluctuations 
shows that all the 'spring constants' (the coefficients of 
a n 2 and b n 2 ) are positive, so the classical paths are local 
minima of the action: 



ggauss ^ 



7T 2 T 

E 

7T 2 T 

Eq 



k 2 + ng\k\ 
k 2 + irg\k\ 



4E 



o Jo 



dr {d T (j>f 



7T5 

4 



sin^(7rT(r - t')) 



ti 2 



^n 2 + ^(\n + k\ + \n-k\-2\k\) 



(4 + b 2 n ). 



(20) 



Thirdly, we believe (although we have not proved) that 
for each value of k, the classical paths </3^'(t) is the unique 
global minimum of Sk ]4>] ■ 

The Fourier modes form a complete basis, so our 
parametrization automatically encompasses the Kor- 
shunov instanton trajectories used by other authors^ In- 
deed, it has been found that instanton techniques, while 
useful in Josephson-junction problems, have limited ap- 
plicability in the current situation because of the failure 
of the 'non-interacting instanton gas approximation'^ 

At first sight the presence of Ttg\k\ in Ea. (|2U|l seems to 
suggest that a large value of g rules out significant finitc- 
winding-number effects. Observe however that the fluc- 
tuation term also depends on k, so we should integrate 
out the residual-phase fluctuations and find whether the 
resulting fc-dependent contribution to the action encour- 
ages or further discourages finite winding numbers. In 
Ea. (|2Ufl . the contribution to residual fluctuations arising 
from the tunneling term vanishes for Matsubara frequen- 
cies n < \k\. As a result, cubic and higher order residual 
phase fluctuations begin playing a role. However despite 



our considering residual fluctuations only to Gaussian 
order, we find that at not too low temperatures, exact 
numerical calculations (discussed in Sec lIV|) are in good 
agreement with our result in Eo. l|23ll and Eq. (|24|l below. 
For very low temperatures, the Gaussian expansion is not 
accurate. 

Let us now integrate out the residual phases, a n and 
b n , from Eo. (|20|l . to obtain the 'effective winding- number 
action' S^ 3 -^ j n the Gaussian approximation. The func- 
tional integration produces a determinant 



D k = \{ 



Tn 2 ? -Kg 
hj i 



2|fc|) 



To eliminate the divergences that are inherent in the def- 
inition of the path integral, we normalize this against the 
determinant Do corresponding to zero winding number. 
The fluctuation contribution to the effective winding- 
number action is thus 



6 



A £ga uss = ^gauss _ ^gauss = ^ ^ 

n 2 T , ^1 ^n 2 + ^ {\n + k\ + \n - k\ - 2\k\) 

= —k 2 + irg\k\ + Y-\n^ 1_^L — (22) 

E 2-2 ^n 2 +«g\n\ 

^ T r ( i + |fc[ + ^E^±M ) r(i + \k\ + X +^±M 



kI + *9\k\ In W1 , i,^ 2 , 3 ( 2 3) 



£ ' r(l + |fc|) 2 T{l + x) 



where x 
leads to 



2irm 



7p. If g is large, Stirling's approximation 



uss _ , " |_2 



fe 2 +7r7 gau SS | fc | 



where 78 auss j g a temperature- and winding-number- 
dependent effective coupling parameter given by 



gauss 

Tk 



<J 



1 



2 ln gE 

irg T\k\ 



const, 



(24) 



where the 'constant' is in O(g ), but may depend on 
k. Observe that ^s auss has a form very similar to the 
fluctuation-dressed coupling constant <? rGn (n), 



9ren{n) 



1 



irg nT 



obtained by disregarding winding numbers, and also that 
to leading order, 7S auss < g lcn . We note in En.lp^l that 
when 7T7 = 1, the phase correlation function evaluated 
ignoring winding numbers, C(r) « g Ien /g = \ - ^ is 
approximately 1/2. Therefore once |r| is large enough 
such that C(t) < 1/2, Ea. (|17|l is no longer accurate, and 
the effect of winding numbers must be considered even if 
j> 1. The criterion 7T7 = 1 sets a crossover temperature 
~ gE^er - * 9 ! 2 marking the onset of tunneling effects. 
This is consistent with numerical calculations that indi- 
cate 



T, = £ /V2(e^ - 1). 



(25) 



For calculating correlation functions at low temper- 
ature, winding numbers as well as higher order resid- 
ual fluctuations need to be considered. In the regime 
7T7 gauss > 1, the contributions from winding number 
trajectories to the correlation function are exponentially 
small, and the correlation function behaves according to 
Ea. (|17fl . In the regime 7T7 gauss < 1 or below and close to 
T*, both even and odd residual fluctuations are impor- 
tant in the action when the winding number is non-zero. 
Due to these nonlinearities, we have not been able to cal- 
culate the correlation function directly beyond quadratic 
order. However we have numerical results (see below) 
for contributions from residual fluctuations around the 
winding number trajectories. 



Finally, we note here that we have studied the model 
in Eq.© from another direction, namely, by doing per- 
turbation theory in powers of g, as opposed to powers of 
l/g as in this paper. The effective action and correlation 
functions thus obtained will be plotted in the results sec- 
tion of this paper for comparison, but the details will be 
described in another paper, since they pertain mainly to 
the small-g limit which is not the focus of this paper. 



IV. NUMERICS 

We now describe our numerical methods. Since the 
phase <p(r) is a real scalar field and the action S[ip] 
is a real functional, the system can be studied using 
path integral Monte Carlo simulation (see, for instance, 
R.efs. ll6ll7l) . 

We first discuss the discretization, which is an in- 
evitable part of Monte Carlo simulation. The path is 
sampled at N values of imaginary time, Tj — je, where 
j = 0, 1, 2, . . . , N — 1. The kinetic term is discretized us- 
ing the 'primitive approximation'^ that is, by assuming 
that the path interpolates linearly between adjacent sam- 
ple points. This is not as crude as it sounds: for the free 
quantum rotor, the primitive action coincides with the 
exact renormalizcd action obtained by integrating out 
all intermediate phases ip T for r ^ je. The dissipative 
term is likewise approximated by a quadrature formula 
based on bilinear interpolation. The double pole in the 
dissipation kernel, a(r) = sin T vT > causes the integrand 
at 'diagonal' grid points i = j to be indeterminate. To 
deal with this, we transfer the quadrature weight sym- 
metrically off the point i = j onto the points i = j ± 1. 
This is equivalent to the scheme of Ref. [l(| 



^(M) = JE^ ~ Vi+iY 



where 



ngs 2 1 



j = 1 or N - 1 : 
j = 2,...,N-2 



2 tfi-Vi 



§a(je) 
a(je). 



(26) 



(27) 



Actually the criterion for a good discretization of the 
action is not how well it approximates the action, but 
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how well it reproduces the correlation functions. It is 
possible to derive a discrete action that, when used in 
PIMC, will produce a correlation function Cj which is 
exactly equal to C(jj) in the Gaussian approximation. 
However, away from the Gaussian limit, this approach 
did not give a significant improvement over the others. 

The time-step, e — (3/N, is restricted by the small- 
est time scale in the problem, which is the lower cutoff 
of the logarithm in -Fb are (t) , that is, 1/gEo. Ideally e 
should be chosen to be a constant multiple of 1/gEo, but 
the thermodynamic integration method described later 
requires the same e to be used for all values of g. Hence, 
we have used a time-step e < 1/8E . The discretiza- 
tion error in C(/3/2) can be estimated by comparing runs 
with different values of e, and can be 1% or more at the 
largest values of g. This error, however, is simply a ten- 
dency to globally overestimate or underestimate C(r), 
corresponding to a small renormalization of the parame- 
ters Eq and g, and does not affect the conclusions of this 
work regarding the asymptotic behavior of C(r). 

To generate the Markov chain, we use the Hamilto- 
nian Monte Carlo (HMC) methodiiSiiS This is a ver- 
sion of the Metropolis algorithm in which new config- 
urations are proposed by evolving old configurations in 
phase space according to Hamiltonian dynamics, with the 
aid of fictitious momenta. The bias introduced during 
inexact Verlet time evolution is compensated exactly by 
the Metropolis rejection step; the scheme can be shown 
to satisfy detailed balance. During each Verlet trajec- 
tory, the configuration evolves ballistically rather than 
diffusively as in the standard Metropolis method, which 
allows a much faster exploration of configuration spaced 
The HMC method is particularly suitable for smooth ac- 
tions such as Eq. (|26() . where there are no constraints 
or collisions to complicate the time evolution. Although 
the long-ranged term in the action contains a double sum, 
this can be treated using FFT methods, so that the com- 
putational cost of each HMC step is 0(N log N) rather 
than 0{N 2 ). 

The HMC method in its basic form still suffers from 
the problem that the Verlet time-step must not exceed 
the period of oscillation of the short-wavelength Fourier 
components of the path, and thus the long-wavelength 
components take many time-steps to go through one os- 
cillation. The solution to this is, preconditioning. One 
can exploit the freedom in the definition of the fictitious 
momenta: instead of giving all beads the same fictitious 
mass, the Fourier modes are given wavelength-dependent 
masses so that all wavelengths oscillate at the same rate. 
If the action had been quadratic in (f>, then the Fourier 
components would represent independent normal modes, 
and one could achieve complete randomization in a single 
trajectory. 

The multilevel Metropolis algorithm 20 is a strong con- 
tender, and has the advantage that winding-number 
changes can be included naturally in the proposal dis- 
tribution. However, at the time of writing we do not 
know an accurate 'level action', so to achieve acceptable 



acceptance ratios, the multilevel method would have to 
be applied to short sections of path, sacrificing the effi- 
ciency of the FFT method. 

The algorithm automatically makes jumps between dif- 
ferent subspaces. However, in order to do so it must pass 
over an energy barrier between the k = and k — 1 sub- 
spaces proportional to 1/e, corresponding to a phase dif- 
ference of 7r between adjacent time-slices. This is N times 
higher than the typical energy difference between a k = 1 
path and a k = path. Hence, winding-number changes 
can be thermodynamically possible but kinctically hin- 
dered, which is an obstacle to achieving ergodicity. One 
approach we tried is the 'clipped-barrier' scheme, which 
operates as follows. During the calculation of the MD 
trajectory, use a modified fictitious force, corresponding 
to an action in which the barriers have been clipped. 
This makes it easier for the time evolution to climb over 
the barriers. Then, use the true action in the Metropolis 
accept/reject step, which restores detailed balance. We 
found that this method produced shorter autocorrelation 
times. 

Nevertheless, a more attractive approach, which elim- 
inates the above problem entirely, is to perform simula- 
tions at fixed winding number, in separate fc-subspaces. 
We compute the correlation functions at fixed winding 
number, Ck(r), for each k: 



1 cos(<^ r -<p T '), 



where the partition function within each subspace is 



Zi- = e 



-s k 



= / [M 



,-Sk 



The full correlator is given by a weighted average 



-AS fc 



C k (r) 



C(r) 



(28) 



where ASk = Sk — Sq is the relative effective winding- 
number action or 'free energy difference'. This can be 
computed by simple importance sampling^ii^ or by 
thermodynamic integration^ We use the latter method, 
as it is more reliable when AS k is large. Define the ther- 
modynamic function ^k = — which can be obtained 
from Cfc(r): 



* fc = ^| dra(r)(l-Q 



?S?k actually diverges as In N, where N is the Trotter num- 
ber (number of time-slices), but A^, = ^ — \&o is finite 
for large N. The error caused by subtracting two large 
numbers is not too serious in practice. We can now ob- 
tain ASk by integration: 



AS k (g) = AS k (0) 



dg' AV k (g') 
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The effort of computing Ck (t) for many intermediate val- 
ues of g is well rewarded, because one can then compute 
^k, ASk and C(r) for each of these g's with no extra 
work. Besides, ASk(g) is the true effective winding- 
number action, which can be compared directly with the 
analytic estimate A5f auss (g), thus providing a valuable 
check of the numerics, and an indication of the range of 
validity of the analytics. 

All the data presented in the next section were com- 
puted using the fixed-winding-number method. The 
computation took the equivalent of 38 Pentium4 2.0GHz 
CPU-days. 1200 runs were performed for 12 values of 
(3 from 2 to 768, 10 values of g from 0.125 to 8, and 
10 values of k from to 9. During each run, the sys- 
tem was equilibrated for 100 HMC steps and the correla- 
tion functions were subsequently averaged over 250,000 
to 1,000,000 steps. The Verlet time-step was taken to 
be 0.1-0.2 of the maximum allowable time-step, and the 
average length of each HMC trajectory was 10-15 time- 
steps. The Metropolis acceptance ratios were typically 
65%-95%. The autocorrelation time was estimated by 
calculating the autocorrelation of the deviation of C(/3/2) 
from its mean value, and was of the order of 2-10 HMC 
steps. 

The Monte Carlo error can be accounted for as follows. 
When (3 is large (768, say) and g is not too large (0.5, 
say), the phases ip(0) and <p(/3/2) are practically uncor- 
rected, so cos^O) — ip((3/2)) is almost symmetrically 
distributed on [—1, 1], with a variance of the order of 1. 
Assume that C(t) falls to a negligible value at r ~ 4, so 
that the interval [0, 768] consists of about 200 indepen- 
dent blocks. Then, translational averaging reduces the 
variance of the C(/3/2) estimator for a single configura- 
tion to about 1/V2TJ0 w 0.07. A run of 1,000,000 HMC 
steps with an autocorrelation time of 8 steps is equiv- 
alent to 125,000 independent samples. Averaging over 
these reduces the variance to 0.07/V125000 « 0.0002. 
Indeed, the data in Figs. El and exhibit MC noise of 
amplitude e" 8 5 w 0.0002. 



V. RESULTS 

Figure ^ shows the relative effective winding-number 
action, ASk(g)- The numerical results are plotted to- 
gether with the Gaussian approximation, Ea. (|23[) . and 
the result from small-g perturbation theory. There is 
evidently a crossover from one regime to the other. It 
is seen that the error in Ea. (|23|l is indeed of the form 
O(g ) x function of k. 

Figure [3 shows the correlator C(r) vs r. It is possible 
to identify at least three distinct regimes in the behavior 
of C(t): 

a. Exponential regime: For small g and small r, the 
'kinetic' term, Eq.J7Jl, dominates, and the system behaves 
like a free quantum rotor. The correlation function de- 
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2 4 6 8 
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Figure 1: AS vs g for f3 = 16, 64, 768 (top to bottom). Re- 
sistor color codes indicate k = 1, 2, . . . , 9. Solid curves: AS 
obtained by integration of calculated from PIMC. Dashed 
curves: AS evaluated considering only Gaussian fluctuations 
about winding number trajectories. Dotted curves: AS from 
perturbation theory about small-g limit, up to 0(g 2 ). At high 
temperatures (top), PIMC results agree with the Gaussian ap- 
proximation in Eg. 12311 . As temperature is lowered (middle), 
non-Gaussian fluctuations become important, and significant 
deviations from the Gaussian approximation are observed. At 
very low temperatures (bottom), Eq. l23H fails completely. 
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Figure 2: ln[C(r)] vs r for j3 = 768 for various values of g. The curves are symmetrical about j3/2. Black curves: PIMC data. 
Brown: C* cxp . Red: C rcnoxp . Green: C log . Blue: C* alg . Below e~ 8 the data is swa mped by Monte Carlo error. For g > 1, 
C(/3/2) evolves logarithmically (green curve) while tT* < 1, T» = E /^/2(e*3 - 1). For rT* > 1, C(r) obeys a power-law 
(blue curve), C(r) = (T /T*) 2 / sin 2 (nTr) , for any finite value of g. 
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cays exponentially (brown lines in Fig.[5J): 

C cxp (r) «exp (-E 0T ). 



(29) 



For slightly larger r, theory and numerics suggest expo- 
nential decay with a reduced effective Eq (red lines): 



C 1 ' 



5 (t) w exp 



E t 



1 + 2g/n 



(30) 



The exponential regime is not very relevant in the context 
of this paper, which is concerned with large g. 

b. Logarithmic regime: As explained in Eci. (|17f) . 
C(r) is logarithmic in the exponentially wide range 
— fp— <C r <C -J^e^/ 2 . Restoring the correct cutoff in 
the logarithm, 



C log (r) = 1 



1 

■Kg 



In 



2e^ gE suittTt 



ttT 



(31) 



c. Algebraic regime: For r > 



Eq.^nU) be- 



comes negative, indicating complete failure of the ap- 
proximation. In fact, the numerically-obtained C(r) be- 
gins falling below C log (r) once C(r) < 1/2, as nonzero 
winding-number trajectories start becoming important. 
However, this downward trend does not continue indefi- 
nitely, but is arrested by an algebraic decay. It turns out 
that for large r, the PIMC calculations agree very well 
with the following prediction of small-g perturbation the- 
ory: 



sin ttTt 



(32) 



where the function p(g) is equal to 

p(g) = - 1) = 8ng + 4ir 2 g 2 + 0(g 3 ). (33) 

In fact, in this regime, the individual Cfc(r) themselves 
are also given by Ea. H32l) . Thus, although nonzero 
winding-number trajectories occur in thermodynamic av- 
erages with a significant probability, a calculation which 
ignores them will fortuitously give the right answer. 

The g = 2 and g = 3 graphs in Fig. [3 clearly show the 
crossover from the logarithmic regime to the algebraic 
regime. 

Our result in Eq. Q32Jl supports previous findings 2 ^ that 
the long time dynamics (C(r) = (T/T*) 2 / sin 2 (vrTr)) of 
the Coulomb gas model, such as Eq.JSJ), is determined 
by the critical nature of the Ohmic dissipation. This is 
also confirmed in conformal field theory calculations of 
the spin correlation function in the Kondo problem. 22 
Further insight comes from Griffiths' theorem 2 ^ that at 
large t, the correlation function C(r) cannot decay faster 
than the interaction. 

Considering a logarithmic law, 

C(r) = l-(l/Trg)ln(gE T), 
for tT» < 1, and a power-law behavior, 

C(t) = (T/T*) 2 / sin 2 (ttTt), 




0.0001 



Figure 3: C(/3/2) vs for g = 0.5, 1, 1.5, 2, 3 (bottom to top). 
Colors indicate values of g. Solid lines are PIMC results. 
Dotted lines are C log (/3/2). Dashed lines are C alg (/3/2). A 
clear crossover from a logarithmic law to a power law can be 
seen as the temperature is decreased. 



for tT* > 1, we obtain the impurity susceptibility, 

Xdc(T) « X dc(0) - (2M 2 /3T,)(r/T,) 2 , T< T*. 

(34) 

Here Xdc(0) is exponentially large in g, of the order of 
Xdc(T*) evaluated using Eci. lfH^I . 

We also present a few more figures. Figure 01 is a log- 
log plot of C(/3/2) vs (3 for various g, clearly showing the 
crossover from the logarithmic regime to the algebraic 
regime. 

Figure 0] is a log-log plot of x{T) vs (3 for various g, 
confirming that at exponentially low temperatures the 
susceptibility saturates at exponentially large values ac- 
cording to Eq. 1341) . Note that for the larger values of 
g, simulations have not been performed at low enough 
T to observe saturation. The errors in x(T) are smaller 
than the errors in C(/3/2) because X; being an integral, 
is dominated by the behavior of C(r) at small t, which 
is less susceptible to MC error. 

Figure [SI is a phase diagram based on the behavior of 
C{(3/2) as a function of g and T. The dashed and dot- 
ted lines represent smooth crossovers rather than phase 
transitions. Only the g > 1 part of the phase dia- 
gram is relevant to the current paper. The crossover 
between logarithmic and algebraic behavior occurs at 
T* = E / yj2(e 1T9 — 1). 



VI. DISCUSSION 

We have studied the paramagnon contribution to the 
spin correlation function C(t) and the static impu- 
rity susceptibility Xdc(T) of a strongly damped [g > 
1) magnetic defect with XY rotational symmetry in 
a metal close to a ferromagnetic quantum critical 
point. Our analysis shows that quantum tunneling 
(droplet 's spin flip) effects are negligible above an ex- 
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Figure 4: X dc vs T for g = 0.125, 0.25, 0.3, 0.5, 1, 1.5, 2, 3, 4, 8 
(bottom to top). The saturation of impurity susceptibility 
below T* and a T 2 deviation from saturation can be seen in 
all but the top three curves, for which T* is extremely small. 
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Figure 5: Phase diagram for behavior of C(/3/2), /3 = 1/T. 
The algebraic region, C(/3/2) = (T/T,) 2 , occurs at a low 
enough temperature for any nonzero value of g. This phase 
is dominated by spin-flip processes. Droplet fluctuations are 
frozen in the logarithmic phase; C(/3/2) decreases slowly as 
1 — (l/ng)\n(gEo/T). As one approaches criticality (i.e. g 
increases), spin- flip processes occur below exponentially low 
temperatures T < T* = E a /y/2(e*9 - 1). 



ponentially small temperature, T* = Eg/ y/2(e ng — 1). 
In this regime, C(r) and Txdc(T) both decrease loga- 
rithmically as T decreases. At very low temperatures 
T <C T*, C(t) is very well described by a power law, 
C(t) = (T/T*) 2 / sin 2 (ttTt), somewhat analogous to ear- 
lier worksSi^ on the Coulomb gas model and equivalent 
Kondo formulations. The impurity susceptibility Xdc(T) 
saturates to an exponentially large but finite value at 
T = and does not show any anomalous divergence. 



Near T = 0, Xdc (T) = X(fc (0)[l - (2/3) (T/T*) 2 ], demon- 
strating the critical nature of the Ohmic dissipation. 
While our results were obtained for an XY defect, we be- 
lieve they should be relevant to defects with full spherical 
symmetry (possibly with additional modifications due to 
the Berry phase). We stress that there is a qualitative 
difference between the cases of IsingSi and XY symmetry 
due to the critical nature of the long-ranged 1 /r 2 interac- 
tion for one-dimensional (meaning one-time dimension) 
ferromagnetic chains^ 

We considered in this paper the action of Eq.JHJ) and 
associated non-perturbative effects associated with spin 
flips. These effects are dominant if, due to the large ferro- 
magnetic polarizability of the host metal, the corrections 
proportional to g (corresponding to the even powers of 
exchange coupling J) are larger than other Kondo cor- 
rections proportional to odd powers of Jn(e F ). As ana- 
lyzed by Larkin and Melniko\*i the condition to neglect 
odd-power Kondo corrections is g > Jn(e F ). Since we 
consider g > 1, it is sufficient to require just Jn(ep) > 
(1 + Fa) 1 ' 2 , while the above condition g > Jn(ep) is 
satisfied trivially for Jn(ep) <C 1. 

The weak coupling case (<? <C 1) is analogous to the 
overcompensated Kondo problem studied by Nozieres 
and Blandin^ and Abrikosov and Migdal22, where the 
number of conduction electron channels N c h coupling to 
the defect is much larger than 2S. Each conduction elec- 
tron channel corresponds to a different orbital quantum 
number. The coupling constant (Jnp) obeys the scaling 
equation 



d(JriF) 
d\nD 



-(Jn F ) 2 + N ch (Jn F ) 3 + c{Jn F ) 4 • ■ • , (35) 



where D is the bandwidth, and c depends on S but not 
on N c h- The first term in Eq. l|35|l depends on the sign of 
the exchange interaction J, and for J > gives rise to 
the conventional antiferromagnetic Kondo effect with a 
Kondo temperature Tr- = D exp(— 1/ Jn F )- The second 
term describes the scattering of a conduction electron 
from the impurity dressed by a number N c h of closed 
electron loops. If the second term is much larger than the 
first, the coupling constant renormalizes towards zero as 
follows: 



J(D) = 



Jo 



N ch (J n F y 
2 



HD /D) 



(36) 



The above scaling is the same as the renormalization 
of g = (n/2)(SJn F ) 2 /(l + F a ) obtained by Larkin and 
MelnikoAii for g <C 1 by summing over parquet diagrams: 



g(T) 



<j 



HEq/t) 



(37) 



A comparison of Ea. l|35|l and Ea. (|37|l shows that the 
number of channels is N c h = 2/(1 + F a ). As the run- 
ning coupling constant decreases, the usual first Kondo 
term in Eq. l|35|l may no longer be disregarded. If J > 
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(antiferromagnetic), J(D) flows towards a stable multi- 
channel Kondo fixed point given by 

J*n F = l/N ch = {l+F a )/2, (38) 

which corresponds to a Kondo temperature 

T K = Eoe- 1 /^ = E e- 2/{1+F "\ (39) 

If J < (ferromagnetic), the coupling constant scales 
to zero. In that case, as in the ferromagnetic Kondo 
problem, the impurity's susceptibility is expected to obey 
a Curie- Weiss law as T — > 0. 

In this paper we have considered g > 1, and perturba- 
tion theory in 1/g also shows that g flows towards zero^S 
as is evident from Ea. (|16fl . Thus, in both the weak- 
coupling (g < 1) and the strong-coupling (g > 1) cases, g 
is renormalized towards zero until very low temperatures 
below Tk — Eo exp[— 2/(1 -I- F a )], when odd-power terms 
of the coupling J become relevant. 

The overall picture emerging from our analysis is that 
the compensation of the droplet 's moment takes place in 
two basic stages. In the first stage, on which we focus 
in the main text, the paramagnon fluctuations, which 
originally enhance a spin S to the large magnetic mo- 
ment M, are gradually stripped off. This quenching of 
the droplet moment takes place below the temperature 
T*, which in our case is given by « Eq exp[— ng/2] = 
E exp{—(irSJriF/2) 2 /(l + F a )]. At the second stage, and 
an even lower temperature Tk, for J > the usual Kondo 
effect compensates the remaining moment leading to a 
multi-channel Kondo fixed point, while for J < 0, the 
impurity spin becomes free. The multi-channel Kondo 
fixed point could be an artifact of our assuming equal 
coupling of the droplet to all the N c h angular momentum 
channels of the conduction electrons. In reality, the cou- 
pling for each channel may be different, and it is possible 
that the impurity spin will get compensated successively 
in different channels. Maebashi, Miyake, and VaxmaSS 
recently analyzed the problem in the Kondo regime and 
concluded from the scaling equations that the coupling 
constants approach a multi-channel Kondo fixed point. 
They also suggested that at unattainably low tempera- 
tures a crossover happens when a single channel wins out. 
In the last stages of preparation and submission of our 
manuscript we became aware of another work 30 which 



arrived to some conclusions similar to ours in the case of 
O(N) spin symmetry in the large- N limit. 

In this paper we considered a dilute system of im- 
purities, so that mutual interaction among impurities 
can be neglected. The size of a magnetic droplet, L, 
is determined by the dispersion relation of the param- 
agnon modes and the proximity to the critical point, as 
L ~ £o/\/l + Fa- As the quantum critical point is ap- 
proached, the size of the droplets grows, and the system 
can be considered dilute only if the density of impurities 
is n imp <§C £o~ 3 (l + ^a) 3 ^ 2 - We also ignore various anoma- 
lous effects which possibly arise in very close proximity 
to the critical pointful and explore only the nearly ferro- 
magnetic Fermi-liquid regime. 

There are several experimental systems, involving im- 
purities with giant magnetic moments in a nearly ferro- 
magnetic host metal, in which it may be possible to study 
magnetic droplet phenomena systematically. 32 Among 
these are iron (Fe) dissolved in various transition metal 
alloys, 33 nickel (Ni) impurities in palladium 34 35 (Pd), 
and cobalt (Co) impurities in a platinum (Pt) hosti^Si 

There are some close connections between the dynam- 
ics of a magnetic moment with XY symmetry and the dy- 
namics of the electromagnetic phase in quantum dots and 
granular metals. The power law behavior of C (t) at long- 
t (or small temperatures) can be associated with inelas- 
tic cotunncling in the literature of mesoscopic physics , 3 - 
Inelastic cotunneling was long understood to be impor- 
tant at low temperature (T < Eq \n.(g~ 1 )) for weak inter- 
grain coupling (g < 1). Our present work shows that 
this is so even when inter-grain coupling is large, the dif- 
ference being that when g > 1, inelastic cotunneling be- 
comes important only at exponentially low temperatures 
T < . The competition of inelastic cotunneling and the 
Coulomb blockade can lead to interesting consequences 
for the transport properties of a granular metali^ 
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